#!/usr/bin/env python
# coding: utf-8

# ### 1 导入必要的包

# In[84]:


import numpy as np
import matplotlib.pyplot as plt 
import copy


# ### 2 创建数据集

# 根据题目创建一个包含3个类簇的数据集
# 
# 蓝色点个数约为20个，红色点个数约为200个，绿色点个数约为30个
# 
# 不同题目使用的数据集不同，方法相同，根据需要运行相应数据集生成代码

# - 第一题

# In[85]:


a = np.random.uniform(-5,-3,20)
b = np.random.uniform(-1,1,20)
cluster1 = np.array([[x, y] for x, y in zip(a,b)])

a = np.random.uniform(-2.5,2.5,200)
b = np.random.uniform(-2.5,2.5,200)
cluster2 = np.array([[x, y] for x, y in zip(a,b)])

a = np.random.uniform(3,5,30)
b = np.random.uniform(-1,1,30)
cluster3 = np.array([[x, y] for x, y in zip(a,b)])


# - 第二题

# In[86]:


# a = np.random.uniform(5,7,20)
# b = np.random.uniform(0.5,2.5,20)
# cluster1 = np.array([[x, y] for x, y in zip(a,b)])

# a = np.random.uniform(-3,3,200)
# b = np.random.uniform(-3,3,200)
# cluster2 = np.array([[x, y] for x, y in zip(a,b)])

# a = np.random.uniform(5,7,30)
# b = np.random.uniform(-2.5,-0.5,30)
# cluster3 = np.array([[x, y] for x, y in zip(a,b)])


# 合并三个簇为一个完整的数据集

# In[87]:


dataset = np.append(np.append(cluster1,cluster2, axis=0),cluster3, axis=0)


# 可视化数据集(未分类)

# In[88]:


for i in dataset:
    plt.scatter(i[0], i[1],c='black', s=10)
plt.show()


# ### 3 随机生成聚类中心

# In[89]:


a = np.random.uniform(min(dataset.T[0]),max(dataset.T[0]),3)
b = np.random.uniform(min(dataset.T[1]),max(dataset.T[1]),3)
center = np.array([[x, y] for x, y in zip(a,b)])


# 聚类中心可视化

# In[90]:


color_map=['blue','red','green']
marker_map=['s','o','v']
for i,cet in enumerate(center):
    plt.scatter(cet[0], cet[1],c=color_map[i],marker=marker_map[i])
plt.xlim(min(dataset.T[0]),max(dataset.T[0]))
plt.ylim(min(dataset.T[1]),max(dataset.T[1]))


# ### 4 生成距离矩阵

# 计算距离矩阵，根据KMeans定义，需逐个计算现有点与三个中心点的欧式距离
# 
# 现有中心矩阵(3,2)，三个点（每个点有各自的横纵坐标）；
# 
# 点矩阵(250,2)，生成的250个点；
# 
# 目标是(250*3)的距离矩阵，250个点分别与3个中心的距离。

# 利用python广播机制，求欧式距离（纯矩阵计算，不使用循环）
# 
# $dist(X,Y) = \sqrt {\sum\limits_{i = 1}^d {|{x_i} - {y_i}{|^2}} }$

# In[91]:


def compute_distances(center, dataset):
    num_center = center.shape[0]
    num_dataset = dataset.shape[0]  # 250

    # (3, 250) 全0
    dists = np.zeros((num_center,num_dataset))  

    # (3, 250) = (3, 250) + (3, 1) 列广播，同行每列都加相同数值
    dists += np.sum(center**2, axis=1).reshape(num_center,1)

    # (3, 250) = (3, 250) + (1, 250) 行广播，同列每行都加相同数值
    dists += np.sum(dataset**2,axis=1).reshape(1,num_dataset)

    # (3, 250) = (3, 250) - 2* (3, 2) * (2, 250) 矩阵点乘
    dists -= 2*np.dot(center, dataset.T)
    dists = np.sqrt(dists)
    
    return dists.T  # 返回250行，3列的距离矩阵，转置后比较符合直觉


# 计算每个点到每个中心的距离矩阵 `dists`

# In[92]:


dists = compute_distances(center, dataset)
dists


# 找到每个点离哪一个中心最近 `nearest_center`

# In[93]:


nearest_center = dists.argmin(axis=1)
nearest_center


# 将每个点的最近聚类中心点作为其分类，插入到第三列作为标签
# 
# 生成预测标签数据集`dataset_predict`

# In[94]:


dataset_predict = np.insert(dataset, 2, values=nearest_center, axis=1)
dataset_predict


# 通过距离矩阵及数据标签，可视化当前状态
# 
# 离哪个中心点最近，即分为那一类
# 
# 计算每个点离自己分类中心的总距离，即为损失函数

# In[95]:


center


# In[96]:


def plot_current_state(dataset_predict):
    plt.figure()
    for label in range(3):
        plt.scatter(dataset_predict[dataset_predict[:,2]==label][:,0],
                    dataset_predict[dataset_predict[:,2]==label][:,1],
                    c='',
                    edgecolors=color_map[label],
                    marker=marker_map[label])
    for i,cet in enumerate(center):
        plt.scatter(cet[0], cet[1],c=color_map[i],marker=marker_map[i], s=50)
    plt.show()

plot_current_state(dataset_predict)
loss = np.sum(dists.min(axis=1))  # 计算每个点离自己分类中心的总距离
print("当前总距离为：",loss)


# ### 5 迭代求出最佳聚类中心

# In[97]:


def fit_kmeans(dataset_predict, dists, epoch):
    break_k = 0
    last_loss = np.sum(dists.min(axis=1))
    for i in range(epoch):
        for label in range(3):
            simple_list = dataset_predict[dataset_predict[:,2]==label]
            x = np.mean(simple_list[:,0])
            y = np.mean(simple_list[:,1])
            center[label] = [x,y]
        plot_current_state(dataset_predict)
        loss = np.sum(dists.min(axis=1))
        print("总迭代轮次：",i+1,"/",epoch,"总距离：",loss)
        
        dists = compute_distances(center, dataset)
        nearest_center = dists.argmin(axis=1)
        dataset_predict = np.insert(dataset, 2, values=nearest_center, axis=1) 
        
        if loss == last_loss:
            break_k += 1
            if break_k == 3:
                break
        else:
            last_loss = loss
            break_k = 0
        


# 通过epoch设置迭代次数，连续3次距离不下降，则自动结束

# In[98]:


epoch = 20
fit_kmeans(dataset_predict, dists, epoch)

